!
!     trashift     shift the image according to a .TRA file
!     Copyright (C) 2004  Filip Hroch, Masaryk University, Brno, CZ
!
! $Id$
!

       Program ProgTraShift
!
! Experimental version 
!
! Shift frames.
! Using the output of the munimatch (with -t switch), eg. 'TRA'
! files with format:
! image1 image2 scale dx1, dx2 cos(fi) sin(fi) sin(fi) cos(fi)

     implicit none

! Parameters:

!    Version:
         Character(len=80),parameter :: version = 'Trashift v. 0.0.0'

! Initial parameters.
!    Output image name:
       Character(len=256) :: namekombi='trashift.fits'

!    Minimum of values
       Real :: minvalue = 0.0

!    Maximum of values
       Real :: maxvalue = 65535.0

!-------

       Character(len=256) :: sb, sbtra
       Character(len=80) :: buf

       Logical :: simple,extend,anyf,exist
       Integer :: n, nalldate, maxalldate
       Integer :: i, j, istat, ostat ,bitpix, naxis, naxes(2) 
       Integer :: BlockSize,pcount,gcount

! image dimensions
       Integer ::  nx = -1, ny = -1 
! image
       Real, allocatable :: KOMBINE(:,:), IMAGE(:,:)
       Real :: bscale, bzero, s
       Real :: oneexp, exptime
       character(len=80) :: onedate
       character(len=80),allocatable :: alldate(:), alld(:)
       Real :: xx,yy,c,r(2),m(2,2),t1(2),t2(2)
       character(len=1000) :: sb1,sb2
       integer :: ix,jy

!------

       write(*,*) trim(version)

! init sum variables
       exptime = 0.0
       nalldate = 0
       maxalldate = 2
       allocate(alldate(maxalldate))

       ! nacteni jmena souboru
       Read(*,'(A)',ERR=666) sb
       istat = 0

       ! input initialisation of the current file
       Call ftopen(25,sb,0,BlockSize,istat)
       Call ftghpr(25,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,istat)
       if( istat /= 0 ) goto 666
       CALL ftgkye(25,'BSCALE',bscale,buf,istat)
       CALL ftgkye(25,'BZERO',bzero,buf,istat)
       !        BSCALE, BZERO not presented
       if( istat /= 0 ) then
          istat = 0
          bscale = 1.0
          bzero = 0.0
       endif
       call ftgkye(25,'EXPTIME',oneexp, buf, istat)
       if( istat /= 0 ) then
          oneexp = -1.0
          istat = 0
       endif
       call ftgkys(25,'DATE-OBS',onedate, buf, istat)
       if( istat /= 0 ) then
          onedate = ''
          istat = 0
       endif
       nx = naxes(1)
       ny = naxes(2)
       Allocate(   IMAGE(nx,ny) ) 
       IMAGE = 0.0
       Allocate( KOMBINE(nx,ny) )
       KOMBINE = 0.0
       call FTG2DE(25,1,minvalue,nx,nx,ny,IMAGE,anyf,istat)
       if( istat /= 0 ) goto 666
       call ftclos(25,istat)
       if( istat /= 0 ) goto 666

       ! image - transformation input
       sbtra = sb(:index(sb,'.f')+index(sb,'.F'))//'TRA'
       open(27,file=sbtra,status='old',iostat=i)
       if( i == 0 ) read(27,*,iostat=i) sb1,sb2,c,r,m,t1,t2
       if( i /= 0 ) then
          write(*,*) 'No such TRA file or invalid format. Skipped.'
       endif
       close(27)
       if( i == 0 ) then

          ! compute sumation
          do i = 1, nx
             do j = 1, ny
                xx = i - t1(1)
                yy = j - t1(2)
                ix = nint(c*(r(1) + xx*m(1,1) + yy*m(1,2)) + t2(1))
                jy = nint(c*(r(2) + xx*m(2,1) + yy*m(2,2)) + t2(2))
                if( 1<=ix .and. ix<=nx .and. 1<=jy .and. jy<=ny ) then
                   KOMBINE(i,j) = IMAGE(ix,jy)
                endif
             enddo
          enddo
       end if

! form the output image
       ostat = 0
       Inquire(file=namekombi,exist=exist)
       if( exist ) then
         open(26,file=namekombi)
         close(26,status='DELETE')
       endif
       Call ftinit(26,namekombi,1,ostat)
       call ftcopy(25,26,1,ostat)
       call ftphis(26,'Shifted image.',ostat)
       if( ostat /= 0 ) goto 666

       Call ftp2de(26,1,nx,nx,ny,KOMBINE,ostat)

       CALL ftclos(26,ostat)
       if( ostat /= 0 ) goto 666

! print output info
       write(*,*) 'Output image: ',trim(namekombi)
       write(*,*) 'Dimension: ',nx,'x',ny
!       write(*,*) 'Bye.'
       stop 'Bye.'

666    continue 

       if( istat +  ostat /= 0 )then
        write(*,*) 'The following error(s) has been occured:'
        do
          Call ftgmsg(buf)
          if( buf == ' ' ) exit
          write(*,*) trim(buf)
        enddo
      endif

      stop 'Crached...'

      End       

